function [Rs, Rm, Cm] = getRC(data, rstepOn, rstepDur, rstepAmp, sampRate, skipCm)
% getRC
%
% Calculates Rs, Rm, Cm from current response to a negative V step.
% Input variables: rstepOn and ~Off are in seconds; 
% rstepAmp is now in millivolts.
% Baseline: 50 ms window before step is taken as baseline.
%
% Editing:
% gs april 2005 -- private version for mapAnalysis2p0
% ---------------------------------------------------------------------------
rstepAmp = rstepAmp/1000;
rstepAmp = abs(rstepAmp);

% stimulus info & indexing
rstepOnIndex = round(rstepOn * sampRate) + 1;   % convert to points
rstepOffIndex = round((rstepOn + rstepDur) * sampRate) + 1;

% baseline
b(1) = rstepOnIndex - round(.05*sampRate); % baseline: preceding 50 ms
b(2) = rstepOnIndex - 1;
baseline = mean(data(b(1):b(2)));
b(3) = length(data); % just used for plotting

% step
s(1) = rstepOnIndex;
s(2) = rstepOnIndex + round((rstepOffIndex - rstepOnIndex)*2/3);
s(3) = rstepOffIndex;

% find peak
[peakY, peakX] = min(data(s(1):s(3))); % assumes pulse is hyperpolarizing; response is downward

% steady state: average over last quarter of pulse
steadystate = mean(data(s(2):s(3))); % stop before pulse end to allow for avg effects

% % tau: exp fit to response decay
% % s(3) = s(1)+40;
% x = (0 : 1 : (length(data(s(1)+peakX:s(3)))-1))';
% y = data(s(1)+peakX:s(3));
% 
% eqn = fittype('a * exp(-x/b) + c');
% opts = fitoptions(eqn);
% set(opts, 'StartPoint', [1 1 1]);
% curve = fit(x,y,eqn,opts);
% a=curve.a;
% tau=curve.b;
% c=curve.c;
% tau2 = tau/sampRate;

% new shorter method
skipFactor = 20;
x = (0 : skipFactor : (length(data(s(1)+peakX:s(3)))-1));
y = (data(s(1)+peakX : skipFactor : s(3)));

% convert x and y to col vectors if necessary
[rw, cl] = size(x);
if rw == 1
    x = x';
end
[rw, cl] = size(y);
if rw == 1
    y = y';
end

% set up fit
eqn = fittype('a * exp(-x/b) + c');
opts = fitoptions(eqn);
set(opts, 'StartPoint', [1 1 1]);
curve = fit(x,y,eqn,opts);
a=curve.a;
tau=curve.b;
c=curve.c;
tau2 = tau/sampRate;

% Calculate Rs
Ipeak = abs(peakY - baseline) * 10^-12;
Rs = (rstepAmp / Ipeak) * 10^-6;

% Calculate Rm
Isteadystate = abs(steadystate - baseline) * 10^-12;
Rm = (rstepAmp / Isteadystate * 10^-6) - Rs;

% Calculate Cm
if ~skipCm
    Req = (Rs * Rm)/(Rs + Rm) * 10^6;
    Cm = (tau2 / Req) * 10^12;
else 
    Cm = NaN;
end

% % % Graphics
% figure
% 
% subplot(2,1,1)
% % plot data from baseline start to 20 msec after pulse ends
% pulseData = data( b(1) : s(3) + round(.02*sampRate));
% L = length(pulseData);
% xdata = [0:L-1];
% xdata = xdata / sampRate;
% plot(xdata, pulseData)
% set(gca, 'XLim', [min(xdata), max(xdata)]);
% % plot(peakX, peakY, 'ro');
% lc = .3;
% % horizontal line at baseline
% line(   'XData', [min(xdata), max(xdata)], ...      
%         'YData', [baseline, baseline], ...
%         'Color', [lc lc lc], ...
%         'LineWidth', 0.5, ...
%         'LineStyle', ':'); % , 'Parent', state.uncaging.analysis.axes
% % horizontal line at steady state
% line(   'XData', [min(xdata), max(xdata)], ...      
%         'YData', [steadystate, steadystate], ...
%         'Color', [lc lc lc], ...
%         'LineWidth', 0.5, ...
%         'LineStyle', ':'); % , 'Parent', state.uncaging.analysis.axes
% subplot(2,1,2)
% plot(curve, x, y, 'b.')
% set(gca, 'XLim', [min(x), max(x)])
% legend('hide')
% 
% % plot(curve)
